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ABSTRACT 

Wc determine an expression for the Type I planet migration torque involving a locally isothermal 
disk, with moderate turbulent viscosity (5 x 10~ 4 < a < 0.05), based on three-dimensional nonlinear 
hydrodynamical simulations. The radial gradients (in a dimensionlcss logarithmic form) of density and 
temperature are assumed to be constant near the planet. We find that the torque is roughly equally 
sensitive to the surface density and temperature radial gradients. Both gradients contribute to inward 
migration when they are negative. Our results indicate that two-dimensional calculations with a 
smoothed planet potential, used to account for the effects of the third dimension, do not accurately 
determine the effects of density and temperature gradients on the three-dimensional torque. The 
results suggest that substantially slowing or stopping planet migration by means of changes in disk 
opacity or shadowing is difficult and appears unlikely for a disk that is locally isothermal. The scalings 
of the torque and torque density with planet mass and gas sound speed follow the expectations of 
linear theory. We also determine an improved formula for the torque density distribution that can be 
used in one-dimensional long-term evolution studies of planets embedded in locally isothermal disks. 
This formula can be also applied in the presence of mildly varying radial gradients and of planets that 
open gaps. We illustrate its use in the case of migrating super-Earths and determine some conditions 
sufficient for survival. 

Subject headings: accretion, accretion disks — hydrodynamics — methods: numerical — planctdisk 
interactions — planets and satellites: formation — protoplanetary disks 



1. INTRODUCTION 

Important interactions occur between a young planet 
and its surrounding gaseous disk. Such interactions are 
generally strong and can lead to planet migration on rel- 
atively short timescales (Goldrcich & Tremaine 1980; Lin 
& Papaloizou 1986; Ward 1997, see also review by Lubow 
& Ida 2010). The shortness of the migration timescales 
may be a problem, since they can be shorter than planet 
formation timescales (> 10 6 years) predicted by the core 
nucleated accretion models of giant planet formation 
(e.g., Lissauer et al. 2009, and references therein). 

There are distinct regimes of planet migration. The 
Type I regime applies to low-mass planets whose tidal 
torques are too weak to open a gap in the disk. In this 
case, the planet is fully embedded in the disk that is 
nearly unperturbed by the presence of the planet. Sev- 
eral calculations have been carried out to determine a 
Type I torque formula for disk-planet interactions, some 
of which are described below. The tidal torque on the 
planet involves opposing contributions from the material 
outside and inside the orbit of the planet. As a result, the 
net torque on the planet depends on both the values of 
the gas properties (such as density and temperature) and 
their radial gradients near the orbit of the planet. Sev- 
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eral studies have investigated the possibility of stopping 
migration by trapping the planet in regions where the 
radial gradients of disk properties may be strong enough 
to counteract the general tendency to inward migration. 
Such disk regions may occur where the density and tem- 
perature vary rapidly in radius due to changes in disk 
opacity or turbulent viscosity (e.g., Menou & Goodman 
2004, hereafter MG04; Matsumura et al. 2007). Disk 
temperature gradients are modified near a planet due to 
enhanced stellar exposure and shadowing effects caused 
by the variations in the vertical thickness of the disk. 
The modified temperature structure due to this effect 
may also act to somewhat slow migration (Jang-Condell 
& Sasselov 2005). 

Most torque calculations to date have assumed that 
the disk is locally isothermal and the disk viscosity is 
not very low (a > 0.001). Under the locally isother- 
mal assumption, the gas temperature is a fixed imposed 
function of distance from the central star that is inde- 
pendent of height above the disk midplane. This ap- 
proximation is formally justified in disk regions where 
the thermal timescale is shorter than the local orbital 
period. This approximation represents a limitation of 
our analysis, but it allows us to directly connect the re- 
sults with those of previous studies that made the same 
approximation. 

Some studies have suggested that important modifi- 
cations to the migration rate can occur if the disk coor- 
bital region responds nearly adiabatically (e.g., Barutcau 
& Masset 2008; Paardekooper & Papaloizou 2008) or 
if the disk turbulent viscosity is very low (Ward 1997; 
Rafikov 2002; Li et al. 2009). In the nearly adiabatic 
case, Type I outward migration appears possible if (1) 
the turbulent viscosity is sufficiently high to allow the 
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corotation torques to be effective, (2) the disk thermal 
timescale near the planet is longer than the time for the 
gas to make U-turns in its horseshoe orbit near the planet 
(of order the planet's orbital period), (3) the disk thermal 
timescale is shorter than the gas libration timescale in the 
coorbital region, and (4) the disk radial entropy gradient 
is negative. The radial entropy gradient is zero if the 
surface density is oc r -3 / 2 , as in the standard Minimum- 
Mass Solar Nebula (Hayashi 1981), and the temperature 
is oc r -1 , as it would be the case in a disk with constant 
opening angle. 

A weak turbulent viscosity regime (a < 10 -4 ) may oc- 
cur in "dead zones" of disks where gas ionization levels 
may be too low to sustain magnetically driven turbulence 
below the disk's outer layers (e.g., Gammic 1996). Un- 
der such conditions, two-dimensional simulations suggest 
that the migration of a ~ 10 Me planet can be halted by 
the strongly nonlinear effects of waves launched at Lind- 
blad resonances that result in shocks (Li et al. 2009; Yu 
ct al. 2010). At sufficiently low levels of turbulent vis- 
cosity (a < 10~ 5 ), vortices can form and collide with the 
planet, producing a jumpy and chaotic form of migration. 

We restrict attention in this paper to locally isother- 
mal disk regions with moderate levels of disk turbulent 
viscosity (5 x 10 -4 < a < 0.05). For such cases, the 
well known disk torque formula of Tanaka et al. (2002, 
hereafter TTW02) describes the situation of a disk that 
is spatially isothermal (the temperature is radially and 
vertically constant) by means of a linear analysis of the 
three-dimensional dynamical equations. Various three- 
dimensional nonlinear simulations have obtained good 
agreement with these migration rates (e.g., D'Angelo 
et al. 2003; Bate et al. 2003). However, there has not 
been a systematic analysis, based on three-dimensional 
simulations, to study the effects of temperature and den- 
sity gradients on migration, which we provide here. 

Others, such as MG04, generalized the three- 
dimensional analytic disk torque calculation of TTW02 
to locally isothermal cases where radial disk temperature 
variations occur. The analysis of MG04 was done using 
two-dimensional resonant torque expressions, together 
with a softened potential to mimic the effects of the 
third dimension (perpendicular to the disk midplane). 
This calculation included only the effects of Lindblad 
torques and neglected the corotation torques. However, 
upon closer inspection of the three-dimensional Lindblad 
torque contribution of TTW02, we find that it disagrees 
with the results of MG04, applied to a fully isothermal 
disk, in important ways. This discrepancy suggests that 
the torque dependence on disk gradients and the trapping 
mechanisms in general require further examination. One 
goal of this paper is to determine the torque formula for a 
planet embedded in a three-dimensional, locally isother- 
mal disk by means of numerical simulations that includes 
effects of density and temperature gradients, as well as 
the coorbital torque (see § 2). This analysis is described 
in § 3. 

The distribution of torques per unit disk mass, as 
a function of radius, provides a description of where 
the disk contributions to the planet torque occur (e.g., 
D'Angclo & Lubow 2008, hereafter DL08). Such 
torque distributions are useful in the application to one- 
dimensional disk evolution models that account for plan- 
ets, including cases where the planets open gaps in the 



disk, the so-called Type II regime of planetary migration 
(see, e.g., Ward 1997, and references therein). Current 
multi-dimensional hydrodynamical codes cannot prac- 
tically follow the disk evolution over long times, such 
as disk viscous times, while one-dimensional codes can. 
The torque density distribution is typically estimated by 
means of the impulse approximation (Lin & Papaloizou 
1986). But that distribution contains a long tail that is 
unphysical, particularly near the star (e.g., Veras & Ar- 
mitage 2004). Another goal of the paper is to provide a 
more accurate formula for the torque density distribution 
that can be applied to such and similar problems, again 
based on three-dimensional simulations. These distribu- 
tions are described in § 4. Some applications of these 
formulae to migrating super-Earths are discussed in § 5. 
Extensions of the torque density formula to situations 
of mildly varying disk gradients and gap-opening plan- 
ets are presented in § 6 and § 7, respectively. Section 8 
contains the summary of the results. 

2. DESCRIPTION OF SIMULATIONS 
2.1. Disk Model 

The simulations are carried out along the lines of DL08. 
We use a spherical polar coordinates system {O; R, 9, </>}, 
whose origin O is located at the star-planet center of mass 
and which rotates at the same angular velocity as the 
planet orbits about the star. The planet's orbit is circular 
and has zero inclination relative to the disk midplane 9 = 
7r/2. Since the disk is assumed to be symmetric relative 
to its midplane, only the disk's northern hemisphere, 9 < 
7r/2, is explicitly simulated. 

We assume that the disk's gas is locally isothermal and 
that the pressure, p, can be written as 

p(R,9,cj>)=p(R,6,<t))(? s (r), (1) 

where p(R, 9, (f) is the mass density. The sound speed of 
the gas, c s (r), is assumed to be a function of cylindri- 
cal radius r = Rsm9. From the equation of state for a 
perfect gas and Equation (1), it follows that the temper- 
ature T oc pj p = c 2 s . Therefore, to determine the effects 
of temperature gradients on disk torques, we vary T in 
radius by varying the gas sound speed, which is taken to 
be a power law in 7'. The hydrostatic disk scale-height is 
defined as 

H(r) = (2) 

where fix is the Keplerian angular velocity. Notice that, 
since the vertical motion and structure of the disk is 
directly computed from the three-dimensional Navier- 
Stokes equations, the actual scale-height of the disk may 
deviate from the quantity H(r). For brevity, however, 
we will refer to H(r) simply as disk scale-height. 

The disk extends in azimuth over 27r radians around 
the star. In models with planet-to-star mass ratio 
M p /M s < 3 x 10~ 4 , the radial disk region extends from 
0.4 a to 2.5 a or from 0.5 a to 1.55 a, where a is the orbital 
radius of the planet. In models with planet-to-star mass 
ratio M p /M s > 3 x 10 -4 , the disk extends from 0.4a 
to 4.0 a in radius. In the ^-direction, the disk domain 
contains from a few to many disk scale-heights above the 
midplane, depending on both r and H(r). 

The initial mass density, p, is independent of the az- 
imuth, (j), and has a Gaussian profile in the meridional 
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direction, 6. In the radial direction, p oc E/_ff and the ini- 
tial surface density, E, is usually taken to be proportional 
to a power of r. In § 6, we consider a more general case. 
The initial flow field has radial and meridional compo- 
nents u/f = —3 v / (2R), where v is the kinematic viscosity 
of the gas, and ug = 0. The initial azimuthal velocity is 
= r [17 — 17k (a)]- The initial (unperturbed) rotation 
rate 17 accounts for the effects of pressure support pro- 
duced by both density and temperature gradients. For a 
pressure supported disk 



n 2 = 17| 



1 

rp 



dp 
dr 



(3) 



Although 17 has a dependence on distance from the mid- 
plane (see, e.g., Equation (4) of TTW02), we adopt as an 
initial condition its value at the midplane that is given 



o 2 = i7| 



(4) 



In the above equation, we set /3 = — dlnE/dlnr and 
£ = — dlnT/dhir. Notice that, in an infinitesimally thin 
(i.e., two-dimensional) disk, Equation (3) involves the 
surface density and a vertically integrated pressure and 
thus Equation (4) becomes 17 2 = 17 K [1 - (/3 + ()(H/r) 2 }. 

We typically impose a kinematic viscosity proportional 
to the inverse of the initial surface density distribution, 
so that the mass accretion rate through the disk (oc v E) 
is nearly constant in radius and the disk can achieve a 
steady state. In such cases, the initial imposed value of 
the parameter (3 is therefore preserved in time. In the cal- 
culations discussed in § 3, v(a) = 10~ 5 a 2 17k (a), which 
corresponds to a turbulent viscosity parameter a = 0.004 
at the orbital radius of the planet when H(a)/a = 0.05. 
However, we also report on models with lower (as low 
as a ~ 5 x 10 -4 ) and higher kinematic viscosity, up to 
p(a) ~ 10~ 4 a 2 Qk(g) ora~ 0.05 (see also § 4). 

We consider several values of the planet-to-star mass 
ratio, M p /M s . For the Type I cases discussed in § 3, 
we mainly concentrate on M p /M s = 3 x 10~ 6 (or 1 Me 
planet if AI S = 1M ). For Type II cases discussed in 
§ 7, we consider planets up to 2 Jupiter masses (M p = 
2 x 10~ 3 M s ). The gravitational potential of the planet, 
$ p , is smoothed over a fraction e of the Hill radius of the 
planet i?H = a [M p / (3M S )] and is given by 
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GM n 



Vs 2 



e 2 R 2 l 



(5) 



where S is the distance from the planet. We generally 
apply a smoothing (or softening) length with parameter 
e = 0.1. In the low- mass regime, another relevant length 
scale is the Bondi radius Rb = GM p /c 2 (e.g., Masset 
et al. 2006; DL08), which can be shorter than Rb- In 
these situations, the smoothing length we apply is of or- 
der i?B or smaller. 

The softening length plays a much weaker role in three 
dimensions than it does in two-dimensional, since the 
vertical thickness of the disk effectively smooths the po- 
tential of the planet. We find that our results are insen- 
sitive to changes in values of the softening length. For 
cases with an M p = 5 x 10 -6 M s planet, changes in the 
softening length by a factor of 10 produced a change in 



the net torque much less than 1%. For cases with an 
M p = 10~ 3 M s planet, changing £ by a factor of 2 re- 
sulted in a net torque variation of about 1%. Further 
details are given in the Appendix. 

2.2. Numerical Method 

The Navier-Stokes equations that describe the motion 
of the disk's gas are written in terms of specific momenta, 
pun, pRug, and pr [u$ + r 17k (a)] (see, e.g., D'Angelo 
et al. 2005), and are solved by means of a finite-difference 
code that computes separately the spatial integration of 
advection and source terms (apparent forces, pressure 
and gravity gradients, and viscous stresses) using an 
operator-splitting scheme (e.g., Zieglcr & Yorke 1997). 
The algorithm employed in the code is second-order ac- 
curate in space and semi-second-order accurate in time. 
Comparisons on disk-planet interaction problems of this 
code with various other codes have been carried out in 
several studies (Kley et al. 2001; de Val-Borro et al. 2006; 
Masset et al. 2006; de Val-Borro et al. 2007). 

The code uses a grid with constant spacing between 
grid points, in each coordinate direction, to discretize 
the momentum equations and features grid-refinement 
capabilities via multiple nested grids (D'Angelo et al. 
2002, 2003). This technique allows us to increase the 
volume resolution by a factor 2 3 for each added level of 
grid nesting. Since we are particularly interested in the 
flow dynamics in a radial region spanning several times 
H on either side of the planet's orbit, to increase the 
numerical resolution we use nested subgrids that extend 
2tt radians in azimuth and cover the whole domain in 
the 0-direction. We employ a variety of grid systems 
that generate a linear resolution AR w a A8 rj a Acf>, in 
the disk region \R— a\ < 4H(a), ranging from 2 x 10 -3 a 
to 5 x 10~ 3 a. In the cases discussed in § 4 that consider 
disk aspect-ratios H/r < 0.01, the grid resolution in the 
region of interest is always such that H(a)/ AR > 10. 

To quantify resolution effects on the radial distribu- 
tions of torques, we perform a convergence study that 
is presented in the Appendix. The results indicate that 
discrepancies are at a level of 1% or better. 

Either nonreflccting (Godon 1996) or wave damping 
(de Val-Borro et al. 2006) boundary conditions are im- 
posed at the inner and outer radii of the computational 
domain. Reflective and symmetry boundary conditions 
are applied, respectively, at the disk surface 6 = 9 m i n 
and at the disk midplane 8 = ir/2 (see also Masset et al. 
2006). 

The calculations discussed in the paper do not involve 
removal of mass from near a low- mass planet. Accreting 
boundary conditions are instead applied for the larger 
planet masses that we consider, M p > 10 _4 M S . Accre- 
tion is performed within a radius ~ 0.1 i?n of the planet. 
The procedure for mass removal is described in some de- 
tail in § 3 of DL08. When accreting boundary condi- 
tions are applied, the removed mass is not added to the 
planet's mass in order to keep it fixed. Except for the 
calculation discussed in § 4.2, the planet is kept on a 
fixed orbit. 

The disk's evolution is typically followed for about 80 
to about 160 orbital periods of the planet. But a number 
of models are evolved for much longer, up to around 2000 
orbits, in order to check for possible torque saturation 
effects. In models that deal with a planet mass M p > 
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10~ 4 M S , the initial density distribution includes a gap 
along the planet's orbit to account for an approximate 
balance between viscous and tidal torques, which reduces 
the relaxation time toward steady state. Torques and 
torque distributions are evaluated when the flow achieves 
a nearly steady state. 

If we denote with dt g the gravitational force exerted 
on the planet by the mass p dV within a grid element of 
volume dV, the associated torque acting on the planet is 
R p x dig, where R p is the position vector of the planet. 
Direct summation of these elemental torques gives the 
total torque. Radial torque distributions are obtained 
by summing elemental torques over the meridional and 
azimuthal directions, as outlined in DL08. Both total 
torques and torque distributions are averaged in time 
over one orbital period of the planet. 

3. TYPE I TORQUE FORMULA 

3.1. Comparison of Three-dimensional and Softened 
Two-dimensional Torques in a Fully Isothermal 
Disk 

Three-dimensional linear analytic calculations of the 
Type I migration rates were carried out by TTW02. 
They provide torque formulae for a disk where the gas 
sound speed is strictly constant, i.e., for an isothermal 
disk. The torque exerted on a planet, including con- 
tributions from Lindblad resonances and the corotation 
resonance, is given by 

T™ = -(1-36 + 0.54/3) 

x S („)nW0i) J (!) 2 , (6, 

where the (azimuthally averaged) gas surface density S, 
the disk rotation rate H, and the disk vertical scale-height 
H are evaluated at the orbital radius of the planet, a. 
We recall that j3 = — dlnS/dlnr. Equation (6) does not 
include the effects of temperature gradients, since the 
disk is assumed to be fully isothermal. 

The corotation resonance contains the complication 
that it can be saturated (have zero strength). In the 
absence of irreversible processes, the flow in the coor- 
bital region follows closed horseshoe and tadpole or- 
bits in the corotating frame of the planet. The steady 
torque associated with such closed orbits is zero (sat- 
urated). Turbulent viscosity introduces an irreversible 
behavior that results in a nonzero torque. The effects 
of turbulent viscosity are stronger for lower mass plan- 
ets (i?H < H), since the radial width of the horseshoe 
orbit region, w « 2a y/ {M p / M s ){a/ H) (see Masset et al. 
2006) , is narrower and the turbulence diffusion timescale 
across this region is therefore shorter (Ward 1992). 

The tables of TTW02 describe the various contribu- 
tions to the torque in the case that both corotation and 
Lindblad resonances contribute in a fully isothermal disk, 
based on linear theory. We infer from these tables that 
the three-dimensional torque contribution due to Lind- 
blad resonances is 

T T L TW = - (2.34 -0.10/3) 

xS( a )^(a)a 4 ^ 2 (|) 2 . (7) 



Mcnou & Goodman (2004) determined the three- 
dimensional torque for a locally isothermal disk (i.e., 
the temperature depends on r) on the basis of two- 
dimensional linear calculations with a softened potential 
(of the type introduced in Equation (5)) to model the 
effects of the third (vertical) dimension. Only Lindblad 
resonances are included. The softening length was cho- 
sen to be H . As they note, it is not clear what value is 
most appropriate within a range of that order. From Fig- 
ure 3 of MG04, we infer that the Lindblad-only torque is 
given by 

Tta = -(0.80 - 0.77/3+ 1.12C) 

xS(a)^(a)a 4 ^ 2 (-|) 2 . (8) 

The subscript emphasizes that this result derives from 
linear calculations that employ a softened gravitational 
potential. We recall that £ = — din T/d In r. As is ap- 
parent, there is substantial disagreement between Equa- 
tions (7) and (8) for £ = 0, where they should agree, in 
particular for the coefficient of /3. Such gradient terms 
can play an important role in describing the possible 
trapping of a planet in regions of density and temper- 
ature radial variations. 

The effects of the density gradient (the /3 term) on the 
torque enter in two ways. With a negative radial gradi- 
ent of the surface density (/3 > 0), the density interior to 
the orbit of the planet is higher than the density exterior 
to it. As a result, the inner Lindblad torques are fa- 
vored over outer ones, giving rise to a outward (positive) 
torque contribution. On the other hand, the negative 
density gradient slows the disk rotation, due to outward 
pressure forces (see Equation (3)). This effect acts in the 
same sense as a drag term would in causing inward mi- 
gration. The degree of cancellation of these two effects is 
critical in determining the net torque contributions by a 
density gradient. In particular, MG04 concluded that the 
cancellation is nearly complete in two dimensions (i.e., 
when an unsoftcned potential is employed), but not in 
three dimensions. However, TTW02 found nearly exact 
cancellation in three dimensions (see Equation (7)). 

There are at least three possible explanations for the 
discrepancy between the Lindblad-only torque of Equa- 
tions (7) and (8) with £ = 0. One possibility is that the 
softened potential applied to a two-dimensional calcula- 
tion, Equation (8), does not adequately represent the full 
effects of three-dimensional torques, Equation (7). The 
analytic calculation of TTW02 suggests that a more com- 
plicated disk response may occur in a three-dimensional 
disk, especially in the presence of gradients, and that 
this response may not be expressible through the use of 
a softened two-dimensional potential. A second possibil- 
ity is that the local two-dimensional torque formula used 
by MG04, from Ward (1997), is not sufficiently accurate 
to determine the effects of gradients on the torque. This 
may be the case since the two-dimensional (no softened 
potential) Lindblad-only torque formula of TTW02 does 
not agree with that of Ward (1997) for a fully isothermal 
disk. The former has a coefficient of (3 equal to « 1.5, 
while the latter has a much smaller value (nearly zero). 
A third possibility is that the three-dimensional Lind- 
blad torque of Equation (7) does not provide an accurate 
description of the Lindblad torque that would occur in 
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the absence of a corotation torque. Equation (7) does 
describe the Lindblad torque contribution in the pres- 
ence of a coorbital torque. The Lindblad and corotation 
resonances can overlap and may affect each other, even 
in the linear theory. However, for /3 = 3/2, the coro- 
tation torque nearly vanishes in the three-dimensional 
calculation and the effects of resonance overlap can be 
largely ignored. (For j3 = 3/2, there is a small, ~ 1%, 
nonvanishing corotation torque contribution caused by 
three-dimensional effects.) For (3 = 3/2, the torque 
predicted by the three-dimensional isothermal analysis, 
Equation (7), has a dimcnsionlcss torque that evaluates 
to —2.2, while the corresponding result for the softened 
two-dimensional potential in an isothermal disk in Equa- 
tion (8) evaluates to +0.36. This large difference in pre- 
dicted torques (and directions) leads us to conclude that 
the resonance overlap, which should be absent in this 
comparison, is unlikely to explain the discrepancy be- 
tween the calculations. 

3.2. Three-dimensional Torque in a Locally Isothermal 

Disk 

In Figure 1, we plot the total torque on the planet as 
a function of density (left panel) and temperature (right 
panel) gradients (/3 and £, respectively) obtained from a 
subset of our simulations (filled circles, see figure caption 
for details) of an M p = 3 x 1(T 6 M s planet (1M E for 
a solar mass star) in a disk of fixed temperature and 
density at the orbital radius of the planet. In general, 
we consider disks with values of j3 in the range from —2 
to 4 and values of Q in the range from —2 to 2. But a 
few models had gradients outside of these ranges, as also 
indicated in Figure 1. The nonlinear simulation results 
we obtain are well fit to the following torque equation: 

T=- (1.36 + 0.62 j8 + 0.43 C) 

xE(a)^( a ) a 4 ^) 2 (|) 2 . (9) 

The dependencies on the planet-to-star mass ratio, 
M p /M s , and on disk scale- height, H , at the planet's or- 
bit (i.e., on the temperature at the orbit of the planet) 
are tested in § 4 for fixed values of /3 and C- We find 
them to agree very well with the scalings in Equation (9) 
over a broad range of variations. 

If we subtract out the corotation contribution in the 
case of an isothermal disk, adopting the corotation torque 
expression given in TTW02, we obtain a Lindblad-only 
torque of 

7£ = - (2.34 - 0.02/3) 

xE(a)^(a)a 4 ^) 2 (^) 2 . (10) 

In this case, we infer a nearly exact cancellation between 
the two torque contributions by the density gradient, as 
described above. Not surprisingly, this result is close to 
the corresponding expression by TTW02, Equation (7). 

The criterion for torque saturation of Ward (1992) is 
that the libration timescalc r^b (over which the specific 
vorticity gradient is removed) is shorter than the viscous 
diffusion timescalc across the coorbital zone r v d ~ w 2 /v 
(over which time the specific vorticity gradient is reestab- 
lished) . The relative velocity of the fluid on the librating 



orbit farthest from the planet's is « |3u>f2/4|, and thus 
riib « 167ra/(3wil), where Q = O(a). The condition 
7iib < t vc j requires a viscosity such that 




On this basis, we expect that the torques in Figure 1 
involve unsaturated corotation torque contributions for 
the viscosity and sound speeds that we impose (a = 0.004 
and H/a = 0.05 at r = a) on a M p = 3 x 10" 6 M s planet, 
since the characteristic a required for saturation is only 
< 10~ 4 . In deriving the saturation condition above, we 
neglect a numerical factor of order unity multiplying the 
right-hand side of Equation (11) that depends on the 
precise value of w. 

If the coorbital torque were partially saturated, it 
would introduce an additional dependence of the over- 
all torque on planet mass. This dependence on planet 
mass would be inconsistent with the quadratic depen- 
dence predicted by linear theory. The torque we obtain 
in Equation (9) has a quadratic dependence on planet 
mass (see § 4), as expected by linear theory. Therefore, 
our expectation that the torque is not saturated is con- 
sistent with its quadratic dependence on planet mass. 
We provide additional evidence for the lack of saturation 
below. 

To test for saturation further, we run models with 
lower and higher viscosity, for a total variation of a of a 
factor 10. If saturation of the corotation resonance oc- 
curred, we would expect to observe a change in the total 
torque equal to 70% of the unsaturated torque, in a disk 
with zero surface density and temperature gradients (see 
Equations (6) and (7)). For such cases, however, we find 
only small differences (generally a few per cent) in total 
torques from those cases based on the standard viscos- 
ity that we apply. The results for somewhat lower and 
much higher viscosity are also plotted in Figure 1 (as 
empty squares and empty diamonds, respectively). The 
small differences in the figure (among different symbols) 
are not due to significant saturation effects, as indicated 
by the time behavior of the total torque illustrated in 
Figure 2, which is nearly constant over time, within a 
margin of 1% or less. 

If the coorbital torque were to saturate, it would gen- 
erally decline (with some oscillations) from its unsatu- 
rated initial value to zero over several libration timescales 
of the coorbital region. Figure 2 plots three cases with 
turbulence parameter a = 0.002, represented as empty 
squares in the left panel of Figure 1 . The torque is nor- 
malized to its average value between 100 and 200 or- 
bits. Since rub ~ 170 orbital periods in these models, 
the torque evolution illustrated in Figure 2 spans a pe- 
riod of 12 libration timescales, or longer. There is no 
evidence of a substantial decline in the total torque and 
therefore no evidence of saturation. 

In § 4, we provide evidence that the torque behavior 
at even lower (a = 5x 10 -4 ) and higher viscosity (up to 
values of a — 0.05) is intrinsic to the torque density dis- 
tribution, whose shape is subject to rather small changes 
with decreasing or increasing viscosity inside this range of 
the turbulence parameter, a. Therefore, we believe that 
our results apply to a disk where the corotation torques 
are unsaturated (i.e., fully effective). 
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Figure 1. Total torque exerted on a low-mass planet by a locally isothermal disk as a function of surface density (left panel) and 
temperature (right panel) gradients, obtained from three-dimensional simulations. We define (3 = — dlnS/dlnr and f = — dlnT/dlnr. 
The quantity S(r) is the azimuthally averaged surface density of the gas. The filled circles represent simulation results of total torques 
normalized to Ef2 2 a 4 (M p /M s ) (a/H) 2 , in which £, Q, and H are evaluated at the planet orbital radius, a. The solid line represents a 
three- parameter, bi-linear fit to data (the leading term in parenthesis in Equation (9)), including data not displayed in this figure. These 
results correspond to a disk kinematic viscosity equivalent to a = 0.004. Results from models with lower (a = 0.002) and higher (a = 0.02) 
viscosity are also plotted as empty squares and diamonds, respectively, for /3 (left panel) and f (right panel) equal to — 1, 0, and 2. 
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Figure 2. Time behavior of the total torque exerted on a low-mass 
planet by a locally isothermal disk, corresponding to the cases with 
low viscosity (a = 0.002) in the left panel of Figure 1. The slope of 
the surface density is indicated by the parameter /3 in the top-right 
corner of each panel. The temperature distribution is such that 
C = 1. The torque is normalized to 71 v , its average value between 
100 and 200 orbits. The libration timescale in these models is 
T lib ~ 170 orbits (see the text). There is nearly no evolution of the 
torque (variations are below 1%) , as is consistent with it being 
unsaturated. The torque remains unsaturated because T vt j < rub 
(see the text for details). 



The torques shown in Figure 1 disagree with the 
predictions, based on two-dimensional models, of 
Paardckooper & Papaloizou (2009) and Paardekoopcr 
et al. (2010). As we argue in § 3.3, torques computed in a 
two dimensional disk are susceptible to smoothing length 



issues and may not always reproduce three-dimensional 
torques. On the other hand, Masset & Casoli (2010) ob- 
tained a level of agreement with Equation (6) that is sim- 
ilar to ours and a similar variation (< 10%) of torques, 
as reported in our Figure 1 , over the range of viscosities 
we consider and for zero disk gradients. 

Equation (9) that we obtain agrees fairly well with 
the unsaturated torque expression in Equation (6) of 
TTW02, with Q set to zero to match the radial isothermal 
assumption in the latter equation. Tanaka et al. (2002) 
provide the dependence of the torque contributions on 
the radial gradients of the surface density, disk thick- 
ness, and pressure in various tables. One might think 
that these tables can be used to infer a dependence of the 
linear torque on the radial temperature gradient. If these 
tables are used in this way, the resulting torque expres- 
sion is in good agreement with Equation (9). However, 
it is not clear if the result is meaningful. The reason is 
as follows. Equation (11) of TTW02 is appropriate for 
a locally isothermal disk with radial temperature varia- 
tions. But since there is a singularity in the term involv- 
ing dSl/dz that they want to avoid, they set the radial 
temperature gradient to zero (so that d^l/dz = 0) in all 
subsequent equations, including the torque derivations. 
In doing so, they not only discard the dtt/dz term but 
also other terms, including a term that contributes to the 
coorbital torque. This coorbital torque term provides a 
contribution that is independent of the vortensity gradi- 
ent (the usual coorbital torque contribution) and is due 
to the radial temperature gradient (see also Casoli & 
Masset 2009). Therefore, the effects of the temperature 
gradient on the linear three-dimensional torque cannot 
be readily deduced from TTW02. 

The coefficient of the temperature gradient that we 
obtain in Equation (9) is significantly smaller than the 
corresponding coefficient in Equation (8), obtained by 
using a two-dimensional smoothed potential. Further- 
more, the contrast between the results is even stronger 
when comparing the ratio of the gradient coefficients to 
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the term that is independent of j3 and £ in the respective 
equations (i.e., 1.12/0.8 versus 0.43/1.36). The slowing 
of migration (or trapping of planets) by means of the 
opacity variations, envisaged in MG04, relies on having 
a sufficiently steeply declining density profile and a suf- 
ficiently shallow temperature profile (see Equation (8)). 
However, the three-dimensional results (Equation (9) or 
(10)) suggest that this trapping mechanism does not oc- 
cur under such conditions. 

Equations (9) and (10) are derived under the re- 
quirement that surface density and temperature gradi- 
ents (or, more precisely, the derivatives dlnS/dlnr and 
dlnT/dhir) are roughly constant over a radial region 
of about 3H interior and exterior to the planet's orbit, 
where the torque is exerted (see § 4). However, in § 4 
we will develop a formalism for computing torques that 
can be applied in the presence of mildly varying radial 
gradients (see § 6) and of gap-opening planets (see § 7). 

3.3. A Two-dimensional Torque Formula for a Locally 
Isothermal Disk 

Two-dimensional (r-(j)) studies of disk-planet interac- 
tions often smooth the gravitational potential of the 
planet, $ p , over a length of order the disk thickness, 
H. This procedure is meant to account for the vertical 
stratification of disk material in proximity of the planet. 
However, as noted by MG04, it is unclear how well this 
approximation works since three-dimensional resonant 
torques depend on the vertical distribution of density in 
the disk (TTW02). Moreover, two-dimensional simula- 
tions of the Type I regime indicate that torques depend 
on the value of the softening length (e.g., Masset 2002; 
Nelson & Papaloizou 2004). 

We perform two-dimensional simulations using settings 
and parameters described in § 2 and applying a planet 
potential with a smoothing length equal to H, the same 
value of the smoothing length ad opted by MG04. The 
potential is then $ p = —GM p /y/S 2 + H 2 , where S is 
the distance from the planet and M p = 3 x 10 -6 M s . Al- 
though the softening length is introduced as a proxy for 
the vertical disk thickness, it obviously affects the hori- 
zontal gradient of the gravitational potential as well. The 
torque density distribution of non-gap-opening planets 
peaks at a radial distance of about H from the planet's 
orbit (see § 4). At those locations, the smoothed gravi- 
tational force, \d$ p /dS\, is s/2/4 of that resulting from 
an unsoftened, point-mass potential. 

The grid resolution in these calculations is Ar = 
aA(f> = 1.5 x 10 -3 . We consider values of the disk gra- 
dients, (5 and £, in the range from —2 to 2. By fitting 
the results from the two-dimensional nonlinear models, 
we obtain a total torque that can be written as 

75d = - (1-17 + 0.15 + O.22C) 

xS(a)^(a)a 4 ^) 2 (|) 2 . (12) 

We check the scaling of Equation (12) with disk thick- 
ness, H, by performing calculations for given values of /3 
and £ and varying the sound speed at r = a. The scal- 
ing with planet mass was already confirmed by previous 
two-dimensional studies (see, e.g., Nelson & Papaloizou 
2004; Masset et al. 2006). 



Clearly, there are important differences between 
our three-dimensional and two-dimensional results with 
smoothed potential, Equations (9) and (12) respectively. 
In fact, although there is not a wide discrepancy in total 
torque for the case of a spatially isothermal disk (Q = 0) 
with constant surface density (/3 = 0), there are factors 
of 4 and 2 difference in the coefficients of parameters (3 
and £, respectively. 

We did not check how sensitive Equation (12) is to 
changes of the softening length and whether it is possible 
to reproduce more closely the three-dimensional formula 
(Equation (9)) with a different a value of this parame- 
ter. However, there are indications that the coefficients 
of f3 and £ in Equation (12) are highly dependent on the 
softening length (S. Li & H. Li 2010, private communi- 
cation) . 

The two-dimensional analytic linear torque of TTW02 
produces a coefficient of j3 equal to 2.83 (see also the two- 
dimensional linear calculations by Korycansky & Pollack 
1993), which is quite different from our value, although 
the torques are nearly identical for (3 = £ = 0. How- 
ever, TTW02 did not apply any smoothing in obtaining 
their result. Therefore, the likely explanation for this 
discrepancy is the smoothing length H that we adopt. 

4. TORQUE DENSITY DISTRIBUTIONS 

Waves carry torques until they are damped and deposit 
angular momentum in the disk. In order to gain some 
understanding of the radial distribution of torques, which 
is not provided by the integrated torque presented in § 3, 
we analyze the azimuthally averaged torque density as 
a function of radius. This function measures the local 
strength of disk torques and should be insensitive to the 
details of wave damping that determine where the torque 
is exerted on the disk. 

The torque distribution per unit disk mass as a func- 
tion of radius, dT(r)/dM , is defined such that 



T 



POO ^rj- 

2tt -j^j(r) £(r) r dr, 
Jo 



dM 



(13) 



where S(r) is the axisymmetric (i.e., azimuthally aver- 
aged) gas surface density. This quantity provides an 
important diagnostic for the nature of the disk-planet 
torques (DL08; Li et al. 2009). The theory of disk reso- 
nances (e.g., Meyer- Vernet & Sicardy 1987; Ward 1997) 
suggests that the torque distribution for a disk with 
smoothly varying properties should approximately follow 



dT . . (r - a 

dM ir) = ^[^r 



(14) 

where J 7 is a dimcnsionless function and H is evaluated 
at the planet's orbital radius, a. (There are some modi- 
fications to this expression due to resonance widths that 
we ignore.) Notice that, by using Equations (13) and (14) 
along with the approximation £(r) ~ 2(a) and perform- 
ing the transformation r — > r/H, one correctly finds from 
the integration that the total torque scales as (a/H) 2 . In 
the past, function T has in effect been taken to be an in- 
verse power law with distance from the planet that is 
modified close to the planet. The expression is based on 
the impulse approximation (Lin & Papaloizou 1986). We 
provide here a more accurate determination of T . 
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Figure 3. Torque density distributions (torque per unit disk 
mass) in a locally isothermal disk, interacting with a low-mass 
planet, for different values of the disk scale-height at the orbital 
radius of the planet. The ratio H(a)/a is indicated in the upper- 
right corner. Surface density and temperature radial gradients are 
taken to be /J = 0.5 and £ = 1 (see § 2.1). The torque is scaled 
by the quantity Q 2 a 2 (M p /M s ) 2 (a/H) 4 , in which fl = Q(a) and 
H = H(a). 
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Figure 4. Torque density distributions in a locally isothermal 
disk, interacting with a low-mass planet, for different values of 
the planet- to-star mass ratio M p /M s , indicated in the upper-right 
corner. Density and temperature radial gradients in the disk are 
the same as those in Figure 3. Again, the torque is scaled by 
D. 2 a 2 (M p /M a ) 2 (a/H) 4 , where U = U(a) and H = H(a). 

In this section, we describe results for testing the uni- 
versality of function J- over a range of disk parame- 
ters by means of three-dimensional simulations. For this 
purpose, we fix surface density and temperature gradi- 
ents and consider the case of a disk with f3 = 0.5 and 
C = 1- This situation corresponds to a steady-state vis- 
cous disk (accretion rate independent of radius) with con- 
stant aspect-ratio, H/r, and constant viscosity parame- 
ter, a. Function T should then depend only on (r—a)/H, 
with these assumed values of (3 and C- 

We test Equation (14) by considering simulations in- 
volving four different values for H/a and M p /M s and lev- 
els of viscosity in Figures 3, 4, and 5, respectively. These 
figures show the results of dT(r) / dM obtained from sim- 
ulations that are scaled such that the curves in the three 
figures should lie on top of each other, if the functional 
form of Equation (14) is correct. The figures demonstrate 
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Figure 5. Torque density distributions in a locally isothermal 
disk, interacting with a low-mass planet, for different values of the 
viscosity parameter a, indicated in the upper-right corner. The 
torque is scaled as in Figures 3 and 4. Surface density and tem- 
perature radial gradients are given by /3 = — dlnE/lnr = 0.5 and 
C = -dlnT/lnr = 1. 

fairly good agreement with Equation (14) for variations 
of H/r. There are ~ 30% variations (measured at the 
curve peaks) in the scaled values of dT(r) /dM (see Fig- 
ure 3), while the unsealed torque densities differ by a 
factor of ~ 20 4 . There is less than a 2% variation in the 
scaled dT{r) / dM for a factor of 20 2 change in torque den- 
sity due to planet mass changes (see Figure 4). There is 
less than 5% change in the scaled dT(r)/dM for a factor 
of 100 change in a (see Figures 5 and 6). These results 
also add support to our claim in § 3 that the coorbital 
disk torques are unsaturated. The reason is that there 
would be a decrease in the magnitude of net torque at 
about the 30% level at higher viscosity, if the coorbital 
torques were to switch from being saturated to being un- 
saturated (see Equations (9) and (10)). In particular, 
Figure 6 indicates that there is no significant change in 
total torque, over many libration timescales, since frac- 
tional variations stay at levels < 0.5%. The slightly lower 
mass planet, M p = 3 x 10~ 6 M s (1 Me), used in § 3 (ver- 
sus M p = 5 x 10 -6 M s adopted here) is even less prone 
to saturation according to the criterion in Equation (11). 

The torque density distributions in, e.g., Figure 4 
would produce a positive total torque, through Equa- 
tion (13), for a sufficiently negative density gradient 
(P > 0) since the multiplication by £ would make the 
positive peak of the profile more positive and the nega- 
tive peak less negative. However, as indicated in Equa- 
tion (14), function J- depends on the radial gradients 
of both density and temperature. For increasing values 
of |/3 1 and the asymmetry between the positive and 
negative peaks becomes larger. As f3 and Q increase, 
the positive peak becomes shallower whereas the nega- 
tive peak deepens. This effect outweighs that produced 
by the surface density. This is illustrated in Figure 7, 
where dl~(r) / dM is shown for various values of f3 (top-left 
panel) and of C (bottom-left panel). In the right panels of 
the figure we plot the cumulative torque, T(r) , defined as 
the integral on the right-hand side of Equation (13) with 
integration limits from to r (and S oc l/r 13 ). As (3 and 
£ increase, the magnitude of the total torque increases. 
The trend can persist for steeper surface densities, up 
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Figure 6. Top: time behavior of the torque density distribution in a locally isothermal disk with turbulence parameter a = 0.005, 0.001, 
and 5 X 10 — 4 , interacting with a low-mass planet. The time in units of orbital periods is indicated in the upper-right (left panels) and 
upper-left corner (right panels). The torque is scaled as in Figures 3, 4, and 5. Disk's surface density and temperature distributions have 
radial gradients given by f) = 0.5 and £ = 1. Bottom: total torque as a function of time, for the same models as in the top panels, 
normalized to the average torque value, 71 v , between 100 and 200 orbits. In the right panel, filled squares refer to case with a = 0.001. As 
in Figure 2 , the lack of significant time evolution of the torque is consistent with coorbital torques being unsaturated. 



to j3 « 6, the largest radial gradient we consider. The 
total torque tends to decrease (i.e., it becomes less neg- 
ative) for positive radial gradients of surface density and 
temperature, but relatively steep functions seem to be re- 
quired to result in a positive torque (see bottom panels 
of Figures 7 and 8). 

4.1. Analytic Approximations of Torque Density 
Distributions 

Analytic expressions of function J- in Equation (14) 
can be found by fitting data. For this purpose, we as- 
sume that J 7 is a parametric function of j3 and £, of inde- 
pendent variable x = (r — a)/H(a), and has the following 
form 



x tanh (pi — p%x), 



(x-ps) /p 6 



(15) 



and pi = Pi(/3, £), with i = 1, . . . , 8, are the parameters 
resulting from the fit. This functional form is not based 
on a physical model of resonant torques, but rather the 
result of searching among simple combinations of func- 
tions whose parameters have direct physical interpreta- 
tions. So, in Equation (15), parameters p\ and p^ are re- 
lated to the amplitude of the positive and negative peaks, 
respectively, while parameters P2 and p§ are the distances 
of the peaks' positions from the planetary orbit. Param- 
eters P3 and pe regulate the width of the positive and 
negative portions of the function. The hyperbolic tan- 
gent affects the slope with which the profile transitions 
from positive to negative and the ratio pijp?, gives the 
intercept of the function with the x-axis. Three exam- 
ples of the fit are displayed in Figure 8, corresponding to 
gradients specified by pairs the (/3, () = (2, 2) in the top 



panels, (0.5, 1) in the middle panels, and (—2, —2) in the 
bottom panels. In the figure, the rescaled torque density 
distribution dT(r)/dM (left panels) and the cumulative 
torque T(r) (right panels) from a simulations (solid dots) 
are compared against the fitting function (solid lines). 

In Table 1, we give parameters pi for a number of rep- 
resentations of function J 7 , corresponding to various pa- 
rameter pairs (/3,(). We find that expression (15) gen- 
erally produces a very good or good fit to the distri- 
butions obtained from simulations. There are various 
ways to quantify how well Equation (15) fits the data, 
one of the most stringent ways is to compare cumulative 
torques generated by the fitting functions against those 
determined from simulations. Three such comparisons 
are plotted in the right panels of Figure 8, which have 
a level of agreement ranging from a few per cent (top 
and middle panels) to ~ 10% (bottom panel) . The ap- 
proximations of function T by means of Equation (15), 
with parameters given in Table 1, produce asymptotic 
cumulative torques that deviate from simulation data by 
typically less than 10%. 

Parameters in Table 1 can be used to generate an ap- 
proximation of function T for a given pair ((3, Q. For 
example, one can proceed by selecting from the table 
four sets of parameters pi corresponding to f3\, /?2, Ci; 
and C2, such that fa < j3 < j3 2 and Ci < C < (2- These 
parameters pi can then be interpolated at (/3, £), and 
the interpolated values can be used in Equation (15) to 
obtain function T . We use bi-linear interpolations to 
generate torque density distributions dT(r)/dM every 
A/3 = AC = 0.2, which were then integrated according 
to Equation (13) to obtain the total torque as a function 
of the surface density and temperature gradients. By 
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Figure 7. Torque density distributions (left panels) and cumulative torques (right panels) in a locally isothermal disk, interacting with a 
low-mass planet, for different values of the radial gradient of the surface density /3 (top panels) and of the radial gradient of temperature f 
(bottom panels). The values of the parameters ft and £ are indicated in the upper right corners of the left panels. Models in the top panels 
have a fixed temperature gradient C = (i.e., the disk is isothermal), whereas in the bottom panels the gradient of E is /3 = 0. Distributions 
dT(r)/dM are normalized to f2 2 a? (M p / M a ) 2 (a/H) 4 ', whereas cumulative torques are normalized to E Q? a 4 (M v /M a ) 2 (a/H) 2 , in which 
E, Q, and H are all evaluated at the orbital radius of the planet, a. 



fitting these data, the coefficients of f) and ( in Equa- 
tion (9) can be recovered with good agreement, within a 
margin of 10%. 

4.2. A Migration Test Case 

The formulae presented above in this section can be 
readily used to describe planet's migration, as indicated 
by the test case discussed here. We set up a three- 
dimensional disk model with a planet whose mass is 
M p = 10- 5 M S (or 3 M E if M s = 1M ). The disk is 
isothermal (f = 0) and has an initial surface density 
corresponding to a parameter /3 = 2. The slope of the 
azimuthally averaged S(r) is maintained throughout the 
disk evolution by imposing a kinematic viscosity v oc , 
so that the disk is in an approximate steady state. 

The planet is allowed to migrate under the action of 
disk torques after 30 orbital periods. We adopt a treat- 
ment to approximately account for the axisymmetric ef- 
fects of disk self-gravity by forcing the planet to rotate 
at (approximately) the same rate as the disk (see Ap- 
pendix C of DL08). The grid rotates about its origin, 
O, at the same angular velocity as the planet (for de- 
tails, see D'Angelo et al. 2005). Therefore, the planet 
moves only radially through the grid because of orbital 
migration. Gravitational forces exerted on the planet 
are computed via direct summation of elemental forces 
and the equations of motion of the planet arc integrated 
by means of a Bulirsch-Stocr algorithm, which uses an 



adaptive timestep control to achieve a relative accuracy 
of 10~ 5 . 

The evolution of the semi-major axis is illustrated in 
the top panel of Figure 9 as a thick solid line. We use co- 
efficients pi in Table 1, for the pair (/?, £) = (2, 0), to ob- 
tain an approximation of function J- (Equation (15)) and 
then calculate the torque density distribution a"7~(r) / dM 
(Equation (14)). The total torque T(a), as a function of 
the planet's semi-major axis a, is evaluated through the 
integral in Equation (13). Conservation of orbital angu- 
lar momentum imposes that da/dt = 2T{a)/{M p a fix), 
which can be written more explicitly as 

dt K{ \mJ\mJ\hJ 

x J\(^^,P,C^Z(r)rdr. (16) 

The integration limits in the above equation are such 
that r 2 — a S> H(a) and a — r\ ^> H(a). 

A numerical solution of Equation (16) is shown as a 
thin solid line in the top panel of Figure 9. The aver- 
age relative difference between the orbital evolution ob- 
tained from the three-dimensional model and that sim- 
ulated with the formulae presented above is at the 10 ~ 3 
level for 10% variations of the planet's semi-major axis 
(see Figure 9, bottom panel). 
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Figure 8. Left panels: torque density distributions in a locally isothermal disk whose surface density and temperature gradients are 
specified by the pair (/3, f ) indicated in the top-right corner of each panel. The dots indicate the three-dimensional simulation results that 
are rescalcd by Q 2 a 2 (M p / M B ) 2 (a/H) 4 . The solid lines represent fits of Equation (15) to the data points, with parameters pi listed in 
Table 1. Right panels: cumulative torques, rescaled by S Q 2 a 4 (M p / M a ) 2 (a/H) 2 , resulting from dT(r)/dM shown in the left panels for 
both simulation data (dots) and fitting functions (solid line). 



5. APPLICATIONS TO EVOLUTION TRACKS OF 
SUPER-EARTHS 

Torque density distributions arc often used to describe 
the tidal interaction between a disk and embedded bod- 
ies. In this section, we describe an application of the 
formulae provided in § 4 to protoplanctary disks. Simi- 
lar applications to other contexts involving astrophysical 
disks may employ an analogous strategy. 

The synthesis of semi-major axis distributions of ex- 
trasolar planetary systems requires solving for the evo- 
lution of the protoplanetary disk and, simultaneously, 
for the gravitational interaction occurring between the 
planet (s) and the disk. Since these calculations involve 



evolution timescalcs equal to typical lifetimes of proto- 
planetary disks (millions of years) and length scales of 
hundreds of AU, they need to rely on one-dimensional 
models, as a consequence of the limitations of current 
computer capabilities. 

The evolution of a thin, axisymmetric and viscous disk 
that interacts with one (or more) embedded planet (s) 
is described by the following equation (e.g., Lin & Pa- 
paloizou 1986) 

where %(r) = -2Trr 3 vZ dtl/dr = 37rr 2 ^£ft (Lynden- 
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Table 1 

Fit parameters for the function T in Equation (15). 



(AC) 



Pi 



in 



p.', 



P4 



P-, 



PC: 



Pi 



(-2,-2 
(0, -2 
(1,-2 
(2, -2 
(4, -2 



0.0501524 
0.0419354 
0.0379796 
0.0315764 
0.0202990 



0.849572 
0.984212 
1.05644 
1.15319 
1.32162 



0.920263 
0.861833 
0.826359 
0.803355 
0.713699 



0.0233445 
0.0293424 
0.0328917 
0.0368841 
0.0464587 



1.20039 
1.03610 
0.945701 
0.760423 
0.557390 



1.15790 
1.18903 
1.21163 
1.35752 
1.48935 



0.145358 
-0.481275 
-0.936489 
-0.989595 

-1.11322 



3.34263 
4.38744 
5.30156 
4.54574 
3.58864 



(-2,-1 
(0,-1 
(1,-1 
(2,-1 
(4,-1 



0.0466666 
0.0382860 
0.0349835 
0.0309523 
0.0228084 



0.863280 
0.997974 
1.05169 
1.11777 
1.34202 



0.966481 
0.910000 
0.902117 
0.889403 
0.789789 



0.0264784 
0.0328630 
0.0365165 
0.0405477 
0.0507717 



1.13846 
1.00086 
0.922390 
0.840493 
0.629326 



1.12637 
1.15666 
1.17324 
1.22073 
1.32866 



0.300767 
-0.171042 
-0.403456 
-0.593193 
-0.897774 



3.29211 
3.58431 
3.31236 
3.16310 
3.07277 



( 



-2,0 
(0,0 
(1,0 
(2,0 
(4, 



0.0421255 
0.0341462 
0.0303427 
0.0287186 
0.0223466 



0.897070 
1.04357 
1.11959 
1.21232 
1.32467 



1.01784 
0.958885 
0.928609 
0.849933 
0.785512 



0.0304875 
0.0369814 
0.0408822 
0.0447468 
0.0542731 



1.02268 
0.933042 
0.861233 
0.802171 
0.646047 



1.14905 
1.13253 
1.15412 
1.16904 
1.22339 



0.428988 

0.0358698 

-0.181823 

-0.466765 

-0.713550 



2.89123 
3.23369 
3.07328 
3.59154 
2.69379 



(-2,1 
(0,1 
(0.5,1 
(1.5,1 
(2,1 
(4,1 



0.0342374 
0.0317793 
0.0297597 
0.0270259 
0.0255936 
0.0178517 



0.956221 
1.06971 
1.09770 
1.22111 
1.27569 
1.39810 



1.11505 
0.969339 
0.938567 
0.891422 
0.867892 
0.912140 



0.0352915 
0.0401954 
0.0421186 
0.0467822 
0.0489470 
0.0595500 



0.814786 
0.931402 
0.902328 
0.832732 
0.791497 
0.630116 



1.25951 
1.04558 
1.03579 
1.06184 
1.06753 
1.16562 



0.393368 

0.221833 

0.0981183 

-0.226587 

-0.358692 

-0.487347 



2.10245 
4.50245 
4.68108 
4.07528 
3.69761 
1.89963 



(-2,2 
(0,2 
(1,2 
(2,2 
(4,2 



0.0347398 
0.0274646 
0.0244322 
0.0211401 
0.0184703 



0.864893 
1.07401 
1.20668 
1.27928 
1.32365 



1.21887 
1.04454 
0.986380 
0.987525 
1.04665 



0.0346082 
0.0435430 
0.0474291 
0.0518687 
0.0666310 



1.00993 
0.924349 
0.866602 
0.785682 
0.671786 



1.03900 
0.975019 
1.01086 
1.05005 
1.01557 



0.823430 
0.538068 
0.323656 
-0.134179 
-0.446362 



3.99184 
5.37855 
4.93965 
3.71919 
1.84684 
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Figure 9. Top panel: a comparison between semi-major axis evo- 
lutions obtained from a three-dimensional model (thick line, see the 
text for details) and from the integration of Equation (16) (thin 
line). Function T uses parameters p; from Table 1. Bottom panel: 
running-time average of Aa/a, defined as (1/i) L (Aa/a)dt' , where 
Aa is the difference and 5 is the mean value of the two curves in 
the top panel. Data are averaged every one (initial) orbital period. 



Bell & Pringle 1974) is the disk viscous torque and — T(r) 
is the tidal torque exerted by the planct(s) at radius r. In 
Equation (17), it is assumed that the disk's rotation rate 
is unperturbed, i.e., Q = £Ik- Since one can write the sec- 



ond term in square brackets as dT/dr = 27rr£ dT/dM, 
Equation (17) becomes 



4^ 



v ^ 9 U r 9 i v r\ 2S 9T t \ 

(18) 

In case of a multiple-planet system, the second term on 
the right-hand side involves a torque density distribu- 
tion for each planet, and thus the summation dT/dM — 
En 9 Tn/dM. 

The left-hand side of Equations (17) and (18) includes 
the term dT, pe /dt, the mass loss flux from the disk. Here 
we assume a simple model of gas removal due to photoc- 
vaporation by the central star that includes effects of UV 
photons 



dt 



3.7 x 10 



-13 



10AU\ 5/2 



( M<; 



VAU 5 



yr 



(19) 



for r > 10 AU and d£ po /<9i = otherwise. Equation (19) 
applies to a disk surrounding a 1 Mq central star that 
emits ionizing photons at a rate of 10 41 s _1 (see, e.g., 
Martin et al. 2007, and references therein). 

Equation (18) is solved by means of an implicit scheme, 
second-order accurate in time, over a grid that extends 
from 0.01 AU to 1850 AU. The grid is composed of 
10000 elements and uses a constant logarithmic spacing. 
Torque density distributions are obtained from Equa- 
tions (14) and (15), at each timestep, according to the 
interpolation procedure outlined in the last paragraph of 
§ 4.1. The parameter /? required to construct function T 
is computed, every timestep, as an average slope 

P = -7Z TT / ~-J-dr, (20) 



(r-2 - ri) 
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in which the integration limits are r2,i = a±3H(a). The 
orbital radius of the planet is computed by numerically 
solving Equation (16). 

The examples presented here are based on two disk 
configurations, both intended to reproduce the initial 
mass distribution of gas in a Minimum-Mass Solar Neb- 
ula. In the first, the initial disk surface density is taken 
from Davis (2005). At 1AU, Ew llOOgciTr 2 at time 
t = 0. Within 4AU, E oc l/\/r, followed by an expo- 
nential decay at larger radii. The disk contains a total 
mass of ~ 0.02 M© inside of about 40 AU. We assume 
that H/r is constant (i.e., the disk temperature T oc 1/r) 
and that the kinematic viscosity is v = aH 2 Q, where a 
is also a constant (hence v oc \Jr). In the second con- 
figuration, S = 1700 g cm" 2 (1 AU/r) 3/2 at time t = 
(Hayashi 1981). An exponential decay is applied beyond 
about 40 AU so that the disk mass within ~ 70 AU is 
~ 0.02 Mq. The aspect-ratio of the disk is H/r oc r 1 / 4 
(and thus T oc 1/y/r). In this case, we assume that 
a oc y/r, hence the kinematic viscosity v oc r 3 / 2 . 

Figure 10 shows the surface density at various times 
(see figure's caption for details). The top panels refer to 
the first disk's configuration, whereas the bottom pan- 
els refer to the second configuration. Each panel corre- 
sponds to a given value of a at 1 AU (but recall that a is 
constant in the first configuration and only varies with ra- 
dius in the second), which produces initial accretion rates 
on the star between 10 -8 Mq yr _1 and 10 -7 Mq yr _1 . 

The bottom-left panel of Figure 10 indicates that a 
density gap starts to form along the planet's orbit at late 
evolution times. This is due to the small viscosity and 
disk thickness at small radii. A condition for gap forma- 
tion is that the symmetrized one-sided torque, i.e., the 
torque exerted by the planet on disk's gas interior and 
exterior of its orbit, exceeds the viscous torque (e.g., Lin 
& Papaloizou 1986). The one-sided torque is of order 
(a/H)\T\. This can be found by using Equation (14) 
and approximating the integral J \J r s(x)\T,(r)rdr in 
Equation (13) as aHT,(a)J^°\J-s(xJ\dx, where x — (r — 
a)/H(a) and J~s(x) = [J-(x) — F(— x)] /2 is a sym- 
metrized form of function T(x). Therefore, the condition 
for gap formation becomes 

^(a)n 2 (a)a*^\^y >%(a), (21) 

with £ = 27r/ °°|J r s'(a;)|(ix. Recalling the definition of 
viscous torque given above, one obtains the condition 

(S) 2 ^(f) 5 ' 

where the factor £ is typically of order unity. In the 
second disk configuration, this condition is satisfied in 
disk regions interior of a few tenths of AU, but it is never 
fulfilled in the first disk configuration. 

The planet evolution tracks illustrated in Figure 11 as- 
sume that a planetary embryo grows at an oligarchic rate, 
dM p /dt oc M p 2/3 , from 0.1 M E to 5 M E over a ~ 10 year 
period. All tracks start at 13 AU and are computed until 
disk gas within tens of AU has almost entirely dissipated 
and da/dt « 0. The left and right panels of the figure 
refer, respectively, to the first and second disk configura- 
tion. 



Observations show that super-Earths, i.e., planets of 
several Earth masses, are also found within a few tenths 
of an AU from their stars, a likely result of orbital mi- 
gration. This occurrence is predicted in the case of the 
second disk configuration, as indicated by the solid-line 
track in the right panel of Figure 1 1 . The first disk config- 
uration would also predict a similar outcome for a some- 
what smaller turbulence parameter (a m 0.001) or for 
different initial conditions of planet's mass and orbital 
radius (see Figure 12). 

The disk configurations introduced above are also used 
to produce the results illustrated in Figure 12. In both 
cases, a turbulent viscosity parameter a = 0.01, at 1 AU, 
is applied. The disk is allowed to evolve without any 
planet until time to- At time to, we assume that a planet 
of mass M p is located at the orbital radius oq. The disk- 
planet system is then allowed to evolve until negligible 
amounts of gas remain in the disk and da/dt drops virtu- 
ally to zero. The orbital radius attained by the planet at 
the end of the calculation is indicated as Ooo , its asymp- 
totic value. Figure 12 shows function of clq for 
cases in which ty = (left panels), 5 x 10 5 years (mid- 
dle panels), and 10 6 years (right panels). Top and bot- 
tom panels refer, respectively, to the first and second 
disk's initial configuration. In each panel, different lines 
styles correspond to models with different planet masses 
(which is taken to be constant): M p — 0.3 Me (short- 
dashed line), 1 Me (long-dashed line), 3A/e (dot-dashed 
line), and 5 Me (solid line). Clearly, initial conditions 
play a key role in the determination of semi-major axis 
distributions of planet populations. 

6. EFFECTS OF GRADIENT VARIATIONS ON 
TORQUE DENSITY DISTRIBUTIONS 

As discussed in § 1, a number of studies have looked 
into the possibility of stopping migration by trapping 
the planet at disk locations where the radial gradients of 
surface density and/or temperature may be strong and 
vary over short radial distances, as it may occur due to 
rapid changes in disk opacity or turbulent viscosity. 

The results developed in § 4 were derived under the 
condition that surface density and temperature slopes, (3 
and £, are nearly constant over a radial region of about 
3H(a) on either side of the planet's orbit, where most of 
the torque is produced. Here, we wish to show that the 
same formulae arc applicable when disk radial gradients 
vary somewhat over this region. This is because torque 
density distributions respond smoothly to mildly varying 
gradients and can be approximated in terms of mean 
slopes. 

We set up a three-dimensional model with a surface 
density E whose azimuthal average is illustrated in the 
left panel of Figure 13. The profile is such that E oc 1/r 2 
at r < 1 and constant at r 3> 1, as also indicated 
by the long-dashed and short-dashed lines, respectively. 
This setup mimics the conditions of the opacity jumps in 
MG04 for a locally isothermal disk. Details of the slope 
P as a function of radius, in the transition region, are 
shown in the inset. This profile is maintained stable in 
time by imposing a kinematic viscosity v oc 1/E. The 
temperature is constant throughout. In the right panel, 
rescaled torque density distributions dT(r) / dM are plot- 
ted when the planet (M p = 10~ 5 M s ) is located at three 
different radial positions, indicated in the legend. 
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Figure 10. Surface density at various times in years (see legends) of a protoplanetary disk that interacts with a low-mass planet (M p = 
5 Me for t > 10 6 years). The top panels refer to a configuration based on the initial E (solid lines) from Davis (2005), constant a, and 
temperature T <x 1/r (see the text for details). The bottom panels refer to a configuration based on the initial surface density from Hayashi 
(1981), a <x yjr, and T <x 1/y/r . From left to right, the turbulence parameter a at 1AU is 0.003, 0.005, and 0.01 for models in the top 
panels, and 0.0035, 0.005, and 0.007 for models in the bottom panels. 




2 4 6 8 0123456 

t (Myr) t (Myr) 



Figure 11. Evolution of the semi-major axis of a planet that grows at an oligarchic rate (dM p /dt <x M p 2 ' 3 ), from 0.1 Me to 5 Me in 
~ 10 6 years, and interacts with a viscously-evolving and photoevaporating protoplanetary disk. The parameter characterizing disk viscosity 
at 1 AU is given in the bottom-right corner of each panel. The left and right panels correspond, respectively, to the first and second disk 
configuration (see the text for details). The insets show details of the semi-major axis evolution over the first 3 X 10 6 (left) and 2 X 10 6 
(right) years. 
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Figure 12. Asymptotic orbital radius, aoo, as a function of the orbital radius at time to (see the text), for planets whose mass is 
M p = 0.3 Me (short-dashed line), 1 Me (long-dashed line), 3 Me (dot-dashed line), and 5 Me (solid line). Times to are given in the 
legends in units of years. In the top panels, the planet interacts with a protoplanctary disk whose initial configuration is based on the 
surface density from Davis (2005), a constant turbulence parameter a = 0.01, and a temper ature T <x 1/r. In the bottom panels, the initial 
disk configuration is based on the surface density from Hayashi (1981), a = 0.01 \frj\ AU, and T <x 1/^/r. Notice that some panels only 
report two mass cases to improve legibility. 



3.0 




Figure 13. Surface density across a disk transition region (left) obtained from a three-dimensional calculation. The long-dashed line 
is the function 1/r 2 . The transition is produced by choosing an appropriate kinematic viscosity u(r). The inset shows details of the 
transition in terms of the slope /3. The distributions dT(r)/dM (right) are scaled by fl 2 a 2 (M p /M s ) 2 (a/H) 4 . The thin solid curves are 
representations of function T (Equation (15)) with parameters pi (Table 1) interpolated at /3 m 0.2 (profile marked by a filled square) and 
~ 1.8 (profile marked by a filled diamond). The third thin solid curve has a slope parameter = 1. All slopes are mean values computed 
with Equation (20). 
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Figure 14. Migration track of an M p = 10 M s planet (solid 
line) that moves across the surface density transition region il- 
lustrated in Figure 13. The track was obtained by numerically 
integrating Equation (16). The surface density is that displayed in 
Figure 13 and the mean slope /3 (to compute function T) is given 
by Equation (20). Also shown are migration tracks through an 
isothermal disk with constant surface density (short dashed line) 
and with £ <x 1/r 2 (long dashed line). 



The torque density distributions in Figure 13 exhibit a 
behavior reminiscent of that shown in the top-left panel 
of Figure 7 where different, but radially constant, slopes 
ft are imposed. This behavior is also seen by com- 
paring the dot-dashed and short-dashed curves in the 
right panel with the thin solid curves marked by a filled 
square and diamond, respectively. These are torque den- 
sities per unit disk mass obtained from Equation (15) 
and parameters pi (from Table 1) interpolated at slopes 
ft ~ 0.2 (filled square) and ft sa 1.8 (filled diamond). 
Both slopes are obtained as mean values according to 
Equation (20). The torque density distribution produced 
when the planet is located at radius r = 1 (long- dashed 
line) has an intermediate behavior and can be approxi- 
mated by using a mean slope ft = 1, again resulting from 
Equation (20). 

Therefore, if a disk has radial gradients that are mildly 
varying over a region of a few times H, one can still 
apply the expressions derived in § 4 after having intro- 
duced appropriate mean slopes. One important aspect 
to appreciate is the small variation of the torque density 
distribution dT(r)/dM as the density gradient changes. 
Even though ft varies by factors of order unity, the peaks 
of dT(r)/dM are only affected at ~ 10% levels (see also 
Figure 7, top-left panel). The results do not show any 
indication of a migration trapping. 

Figure 14 shows, as a solid line, the migration track of 
an M p — 10 -5 M s planet across the transition region of 
Figure 13. The migration track is obtained by integrating 
Equation (16) (E is that illustrated in Figure 13) and 
using Equation (20) to evaluate the mean slope ft that 
characterizes function J 7 . For comparison purposes, the 
figure also shows the tracks that the same planet would 
follow if £ was constant (short dashed line) and £ oc 1/r 2 
(long dashed line). The solid and_short dashed lines are 
similar until the density gradient ft becomes greater than 
zero. Incidentally, the slope of the long dashed line is 
similar to the other two initially since the reduced density 
value compensates for the larger density gradient (see 



Equation (9)). 

The rapid decay at small radii of the planet's semi- 
major axis seen in the solid and long dashed curves of 
Figure 14 is what one would expect from, e.g., Equa- 
tion (9) applied to a stationary disk: |T(a)| oc a^~@\ 
hence 

,(C-0+i/2). (23) 



da 
~db 



The solid curve is therefore expected to steepen (\da/dt\ 
diverges) as the semi-major axis decreases, whereas the 
short dashed curve is expected to flatten (\da/dt\ tends 
to zero) as a decreases. 

Density transitions occurring at the edges of a "dead 
zone" may be quite abrupt and the slope ft may change 
by large amounts over radial distances of order H (e.g., 
Matsumura et al. 2007). The analysis provided here may 
not apply to those sharp transitions, whose effects on 
torques would require a more detailed analysis. 

7. TORQUE DENSITY DISTRIBUTIONS FOR 
GAP-OPENING PLANETS 

The torque density distributions given in § 4 were de- 
rived in a regime in which the mass density distribu- 
tion along the planet's orbit remains basically unper- 
turbed. In general, this condition ceases to hold when 
tidal torques exceed viscous torques (see Equation (22)), 
which is roughly expressed by the inequality (see § 5) 



( A_\ /2 (Ml 

\3na) \M S 



5/2 



> 1. 



(24) 



For standard disk parameters H/a ~ 0.05, a of a few 
times 10 -3 , and of order unity, Equation (24) re- 
quires a planet-to-star mass ratio M p /M s > 10~ 4 . Gap 
formation can induce average variations in the surface 
density, relative to the unperturbed state, of 2 orders of 
magnitudes when the planet mass grows from one Earth 
mass to one Jupiter mass. However, the corresponding 
variation in the strength of the peaks of the rescaled 
dT(r)/dM is only a factor between 2 and 3, as shown 
in Figure 15. This figure displays torque density distri- 
butions produced by planets of various masses, for cases 
in which i?u < H (top panel) and i?H > H (bottom 
panel). Notice that the radial positions of the peaks are 
located at r — a as ±i?H, when Rn > H (see also DL08). 

The thin solid curve with filled circles in the top panel 
of Figure 15 is function T (Equation (15)) of variable 
x = (r— a)/ H(a) and parameters pi given in Table 1, for 
disk gradients corresponding to ft = 0.5 and £ = 1. The 
same function, but divided by 2 and of variable x = (r — 
a) /Rb, is plotted in the bottom panel. The peak strength 
of the rescaled dT(r)/dM shows only small changes for 
i?H > H. (For i?H ^ H, the torque density scalings can 
be different from those in Equation (14)). The factor 2 is 
given by the left-hand side of Equation (24) for i?n ~ H 
and a more accurate determination of s/^ (see definition 
in § 5). 

For a gap-opening planet, the details of the torque den- 
sity distribution well within the gap region (|r— a\ < Ru) 
arc not very important for computing torques in the 
framework of one-dimensional disk models (of the type 
discussed in § 5) since those portions of dT(r)jdM are 
multiplied by a surface density that is relatively small 
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Figure 15. Torque density distributions in a locally isothermal 
disk for various values of the planet-to-star mass ratio M P /M B , 
indicated in the legends (see the text for further details). The 
disk has parameters H/a = 0.05 and a = 0.004. In both panels, 
the torque is scaled by fi 2 a 2 (M p /M s ) 2 (a/H) 4 . The distributions 
dT(r)/dM are shown for cases where ijjj < H (top) and > H 
(bottom). The short-dashed curve in the two panels refers to the 
same model. See the text for further details. 



(because of mass depletion). For reference purposes, in 
the bottom panel of Figure 15 we also plot as a thin solid 
curve with empty circles the torque density per unit disk 
mass that results from the impulse approximation, as 
implemented by Veras & Armitage (2004). 

8. SUMMARY AND DISCUSSION 

We obtained an expression, Equation (9), for the disk- 
planet torque in three dimensions for a locally isothermal 
disk that contains radial density and temperature vari- 
ations, based on numerical simulations. This expression 
was derived for density and temperature gradients that 
are fairly constant over the torque producing region of 
the disk, which extends roughly over a radial distance of 
3H inside and outside the orbit of the planet, (see, e.g., 
Figures 3, 4, and 5). In the case of a fully isothermal disk, 
the torque expression agrees well with previous analytic, 
linear results obtained by TTW02 (Equation (6)). The 
analysis presented by TTW02 cannot be used to deduce 
the dependence of the torque on the temperature gradi- 
ent (see discussion in § 3.2), which we provided in Equa- 
tion (9). We tested the torque expression for a range 
of disk densities, temperatures, viscosities, and planet 
masses and found good agreement. The results agree well 



with the expected scalings of linear theory with planet 
mass and gas sound speed (see Figures 3 and 4 for which 
linear theory predicts the curves overlap in each figure). 
Although these torque scalings agree with linear theory, 
we cannot claim that Equation (9) with a nonzero tem- 
perature gradient agrees with linear theory because, as 
noted above, there is no expression available for the lin- 
ear torque in a three-dimensional locally isothermal disk. 
These results suggest that substantially slowing or halt- 
ing planet migration by means of smoothly varying but 
large negative density gradients and small temperature 
gradients in a locally isothermal disk is difficult. In ad- 
dition, considerably reducing or stopping migration in 
a locally isothermal disk by means of opacity variations 
(e.g., MG04) or shadowing effects (e.g., Jang-Condell & 
Sassclov 2005) appears unlikely (see § 6). 

We also obtained an expression for the torque den- 
sity per unit disk mass (Equations (14) and (15)). We 
tested the range of applicability of this expression for 
a range of disk temperatures, planet masses, and vis- 
cosities and found good agreement with the expected 
scaling predicted by linear theory. This expression can 
be applied to one-dimensional simulations of disk-planet 
interactions, such as, but not limited to, the cases of 
super-Earth evolution discussed in § 5. The torque den- 
sity expression commonly used is based on the impulse 
approximation and contains a long power-law tail that 
extends far from the planet. Ad hoc modifications have 
been made to make it more physically reasonable (e.g., 
Veras & Armitage 2004). The expression we obtained is 
more accurate and overcomes that difficulty The formal- 
ism presented in § 4 can also be employed for comput- 
ing torques in the presence of mildly varying disk gra- 
dients (see Figure 13) and of gap-opening planets (see 
Figure 15). Moreover, applications to other contexts 
involving tidal interactions between astrophysical disks 
and embedded objects may also be possible. 

The Lindblad-only torque formula of MG04, Equa- 
tion (8), based on the torque density of Ward (1997) 
modified for a planet potential smoothed over a distance 
H, results in a positive torque (outward migration) for 
a steep enough and radially decreasing surface density. 
We notice that, in the fully isothermal case, there is dis- 
agreement between the two-dimensional (no softened po- 
tential) formula of TTW02 (for Lindblad torques only) 
and the results of Ward (1997). At the end of § 3.1, 
we discussed some possible reasons for the discrepancy. 
But, we cannot be certain at this point what the reason 
is. In addition, we found that corotation torques play a 
significant role and, therefore, the use of a Lindblad-only 
torque formula may not be appropriate in the regimes 
that we investigated. 

For zero radial density and temperature gradients in 
the disk, the torques resulting from nonlinear calcula- 
tions of a two-dimensional disk, with a planet potential 
smoothed over a distance H, are in fair agreement with 
the three-dimensional results. However, there is substan- 
tial quantitative disagreement in the effects of these gra- 
dients on the torques (compare Equations (9) and (12)). 
One possible source of this disagreement may be due to 
the effects of the softened two-dimensional potential in 
modifying the torques in the presence of disk gradients. 

The analysis presented here is limited by the range 
of the parameter space we considered, involving a tur- 
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bulcnt viscosity 5 x 10~ 4 < a < 0.05. In addition, it 
is limited by the assumption that the disk thermody- 
namics behaves in a locally isothermal manner. The re- 
sults may not hold under conditions very different from 
those we assumed. For example, if the gas behaves some- 
what adiabatically, the torque can be modified and may 
even result in outward migration (e.g., Baruteau & Mas- 
set 2008; Paardekooper & Papaloizou 2008; Kley et al. 
2009). Under such conditions, nonlinear effects may be 
more important than they are in the locally isothermal 
case. 

In weakly turbulent disk regions (a < 10 -4 ), such 
as in "dead zones", we expect the disk response to a 
planet to be dominated by the highly nonlinear effects 
of waves launched at Lindblad resonances resulting in 
a complex behavior, including the formation of vortices 
(Li et al. 2009; Yu et al. 2010). Such effects produce 
strong changes in the underlying gas distribution near 
a low-mass (~ 5 Me) planet that in turn affect its mi- 
gration rate. Type I planet migration can cease and a 
vortex dominated, jumpy form of migration can prevail. 
The linear torque formulae presented above are inappli- 
cable in such a regime. For a > 5 x 10 -4 , the viscosity 
regime we considered, the nonlinear feedback effects of 
shocks on migration subside due to turbulent diffusion 
that washes out these sharp features in the gas density 
(Li et al. 2009). However, these simulations were per- 
formed in two dimensions and it would be useful to un- 
derstand how three-dimensional effects could alter this 
form of migration. In particular, the formation of vor- 
tices at gap edges may operate differently in three di- 
mensions. 

In highly turbulent disk regions (a > 0.1), there may 
be other types of changes in the properties of the flow 
that can potentially affect disk-planet interactions. Pre- 
liminary results suggest that the shape of the torque den- 
sity distribution may also be altered. Therefore, in such 
disk regions, the torque acting on low mass planets may 
not be represented by the formulae presented here. 

Paardekooper & Papaloizou (2009) performed linear 
and nonlinear two-dimensional calculations of torques 
acting on low-mass planets in fully isothermal disks and 
found that nonlinearities can develop, depending on the 
level of viscosity. For a M p ~ 10 -5 M s planet embedded 
in a disk with H / a = 0.05, they found agreement between 
linear and nonlinear calculations only at very high vis- 
cosity, beyond a ~ 0.1. Those two-dimensional results, 



however, are highly sensitive to the value of the smooth- 
ing length applied in the planet's potential. As discussed 
above, in the fully isothermal case, our three-dimensional 
torques are in good agreement, including the effects of 
density gradients, with the linear theory of TTW02 at 
much smaller viscosity levels, a ~ 0.001. In the more 
general locally isothermal case, we found that the over- 
all scalings of the torque with planet mass and gas sound 
speed agree with the predictions of linear theory. Nonlin- 
ear feedback effects due to partial saturation of the coro- 
tation torque would in general give rise to different scal- 
ings with planet mass (e.g., Ogilvic & Lubow 2003). Such 
effects are not expected to be present in the parameter 
space we investigated (see Equation (11)), as we also veri- 
fied in simulations (see Figures 2 and 6). Other nonlinear 
effects arc possible, especially in other regimes as men- 
tioned above, but we did not see any evidence for them 
to have a significant impact in these three-dimensional 
simulations. 

The analysis provided here may not be applicable to 
cases where substantially large variations of density and 
temperature gradients occur within the disk zone where 
torques are produced, i.e., a few times H or less in radius. 
Such situations may arise near the edges of "dead zones" 
(e.g., Matsumura et al. 2007), although there are limits 
on the slope of these edges imposed by disk instabilities 
(e.g., Yang & Mcnou 2010). It would be desirable to 
explore such situations along the lines of the approach 
taken in this paper. 
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APPENDIX 

A. SENSITIVITY OF RESULTS TO IMPOSED PARAMETERS 

In numerical simulations of the kind discussed in this paper, there is often much emphasis placed on the form of the 
gravitational potential of the planet. In a context in which disk self-gravity is neglected and the "effective" radius of 
the planet is much smaller than the smallest length scale resolved numerically, this potential is that of a dimensionless 
mass M p , i.e., $ p = —G M p / S, where S is the distance from said mass. Since this potential is singular at the position of 
the planet and the planet lies on the computational grid, the singularity could generate numerical instabilities because 
the source terms of the momenta equations could involve the calculation of a diverging gravitational force. (Notice 
that there is no problem with employing a point-mass potential for the gravity field generated by the star, when it is 
located outside of the computational grid.) Therefore, the singularity is usually removed by introducing a modified 
potential as in Equation (5). In this case, the potential is softened over a given length, which is cither a fraction of 
the planet's Hill radius, i?H, or on the order of the local (hydrostatic) disk scale- height, H(a). The second option is 
meaningful only in a two-dimensional (r-0) geometry, where the smoothing length is used, in effect, as a proxy for the 
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Figure Al. Top: torque density distributions in a locally isothermal disk, interacting with an M p 

M p = 10 -3 M s planet (right), for different values of the softening parameter e (see Equation (5)), as indicated in the legends. Torques 
are rescaled by f2 2 a 2 (M p /M s ) 2 (a/H) 4 , where f2 = Q(o) and H = H(a). The curves overlap each other and total torques differ by ~ 1% 
(or less for the lower mass planet). The differences of dT(r)/dM are shown beneath the torque density distributions. On the left, the 
differences are between cases which use e = 0.05 and e = 0.1 (long-dashed line) and cases which use e = 0.01 and e = 0.05 (short-dashed 
line). Bottom: convergence study of dT(r)/dM for the models in the top panels, for cases that use average linear resolutions of 2 X 10 -3 
and 10 -3 (left) and of 4 X 10 -3 and 2 X 10 -3 (right). The agreement is again at the ~ 1% level or better. 



third (vertical) dimension of the disk. We choose the first option and adopt a softening length that is a tenth of i?H • 
The deviation of the softened potential in Equation (5) from the point-mass potential, normalized to the point-mass 
potential, is sa (e 2 /2)(i? H /5 1 ) 2 , which is < 12% for S/Rn > 2s and < 2% for S/Rn > be. We check that our results 
are largely insensitive to the choice of e = 0.1 by performing calculations for a M p = 5 x 10 -6 M s planet in a disk, 
whose viscosity parameter is a = 0.005, and using softening lengths (eRh) with e = 0.05 and 0.01. The resulting net 
torques agree at about the 10 -3 level, implying that they are basically converged with respect to the parameter e. 
We perform a similar experiment with an M p = 10~ 3 M s planet, using e = 0.05, and found discrepancies at the 10~ 2 
level. Torque density distributions are plotted in Figure Al {top panels) for the M p = 5 x 10 -6 M s planet {left) and 
the M p = 10~ 3 M s planet {right), with values of e indicated in the legends. Given the level of agreement, profiles 
corresponding to different softening length cannot be distinguished from one another, or barely so. The differences of 
dT{r)/dM , taking cases with closest values of e (see figure's caption for details), are shown beneath the torque density 
distributions. 

The outcome illustrated in the top panels of Figure Al can be explained with the property of the relative distance 
between the planet and the locations where momenta and gravitational potential are computed. The planet position 
is fixed on the (rotating) grid at a "corner" of a volume clement of the grid (where eight such elements are in contact). 
Momenta are computed at the geometrical center of each of the six faces of a volume element, whereas the gravitational 
potential is computed at the geometrical center of the volume element. Therefore, there is always a buffer distance 
of V3L/2 (neglecting curvature), where L is the average grid spacing, which prevents $ p from diverging, even if 
£ = 0. This argument also applies when the planet migrates since, by virtue of the time-varying rotation rate of the 
reference frame (see § 4.2), it only moves along the radius, R, but not along the azimuth, <f>. Hence, there still exists 
a minimum buffer distance of y2L/2 between the planet position and the center of a volume element. (Notice that, 
on a two-dimensional grid, these buffer distances become shorter and equal to \/2L/2 and L/2, respectively.) 

Numerical resolution is another potential issue that may have an impact on simulation results. In the bottom 
panels of Figure Al , we test for convergence of torque density distributions with respect to numerical resolution. As 
mentioned in § 2.2, the average linear resolution in the disk region of interest {\R — a\ < AH) applied in this study 
ranges from 2 x 10~ 3 a to 5 x 10 -3 a. These grid systems contain typically from 10 to over 40 millions of grid elements. 
The torque density distributions shown in Figure Al {bottom panels) are calculated using average linear resolutions of 
2 x 10 -3 and 10 {left) and of 4 x 10~ 3 and 2 x 10~ 3 {right). The higher resolution cases in these plots are based on 
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grid systems that have nearly 110 (left) and 90 (right) million grid elements. The profiles lie on top of each other and, 
again, can be hardly distinguished. The differences are displayed beneath the plots of dT(r)/dM , showing agreement 
at a ~ 10 -2 level or better. 
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